!
! wcswrite - write well-formated FITS header
!
!
! Copyright © 2011 F.Hroch (hroch@physics.muni.cz)
!
! This file is part of Munipack.
!
! Munipack is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
! 
! Munipack is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with Munipack.  If not, see <http://www.gnu.org/licenses/>.
!

program wcswrite

  use astrofits

  implicit none

  integer, parameter :: ndim = 2
  integer, parameter :: dbl = selected_real_kind(15)
  integer, parameter :: clen = 80
  integer :: ncom, n
  character(len=32000) :: line, record, file='',output = ''
  character(len=clen), dimension(ndim) :: ctype
  real(dbl), dimension(ndim) :: crval, dcrval, crpix, dcrpix
  real(dbl) :: scale, dscale, angle, dangle
  character(len=clen), dimension(:), allocatable :: com
  character(len=clen) :: type
  integer :: istat, eq, cc
  logical :: fcrpix = .false.


  ncom = 0
  n = 0
  ! read input parameters
  do
     read(*,'(a)',end=99,iostat=istat) record
     if( istat /= 0 ) then
        write(*,*) 'Unexpected string on input.'
        stop 21
     end if

     cc = index(record,'#')
     if( cc == 0 ) then
        line = record
     else
        line = record(1:cc-1)
     end if

     eq = index(line,'=')

    ! astrometrical parameters

     if( index(line,'PROJECTION') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) type
        if( istat /= 0 ) then
           write(*,*) 'PROJECTION string read failed.'
           stop 21
        end if
        ctype = type
        

     else if( index(line,'CRPIX') == 1 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) crpix
        if( istat /= 0 ) then
           write(*,*) 'CRPIX parameter read failed.'
           stop 21
        end if
        fcrpix = .true.

     else if( index(line,'DCRPIX') == 1 .and. eq > 0 .and. fcrpix ) then

        read(line(eq+1:),*,iostat=istat) dcrpix
        if( istat /= 0 ) then
           write(*,*) 'DCRPIX parameter read failed.'
           stop 21
        end if

     else if( index(line,'CRVAL') == 1 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) crval
        if( istat /= 0 ) then
           write(*,*) 'CRVAL parameter read failed.'
           stop 21
        end if

     else if( index(line,'DCRVAL') == 1 .and. eq > 0 .and. fcrpix ) then

        read(line(eq+1:),*,iostat=istat) dcrval
        if( istat /= 0 ) then
           write(*,*) 'DCRVAL parameter read failed.'
           stop 21
        end if

     else if( index(line,'SCALE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) scale
        if( istat /= 0 ) then
           write(*,*) 'SCALE parameter read failed.'
           stop 21
        end if

     else if( index(line,'DSCALE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) dscale
        if( istat /= 0 ) then
           write(*,*) 'DSCALE parameter read failed.'
           stop 21
        end if

     else if( index(line,'ANGLE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) angle
        if( istat /= 0 ) then
           write(*,*) 'ANGLE parameter read failed.'
           stop 21
        end if

     else if( index(line,'DANGLE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) dangle
        if( istat /= 0 ) then
           write(*,*) 'DANGLE parameter read failed.'
           stop 21
        end if

     else if( index(line,'FILE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) file,output
        if( istat /= 0 ) then
           write(*,*) 'FILE and OUTPUT parameter read failed.'
           stop 21
        end if

     else if( index(line,'NCOM') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) ncom
        if( istat /= 0 ) then
           write(*,*) 'NCOM parameter read failed.'
           stop 21
        end if
        allocate(com(ncom))
        
     else

        if( ncom > 0 ) then
           n = n + 1
           if( n <= ncom ) com(n) = line(:clen)
        end if

     end if


  end do
99 continue

  if( allocated(com) ) then
     call wcsupdate(file,output,type,crval,dcrval,scale,dscale,angle,dangle,fcrpix,crpix,com)
  else
     call wcsupdate(file,output,type,crval,dcrval,scale,dscale,angle,dangle,fcrpix,crpix)
  end if
  
  if( allocated(com) ) deallocate(com)

  stop 0

end program wcswrite
